{
 "metadata": {
  "name": "",
  "signature": "sha256:568685ce4f1efd25c1c5f2579f75f35270c5a882a9d1ca1fedb550611442d63b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Sparse covariance estimation for Gaussian variables\n",
      "\n",
      "A derivative work by Judson Wilson, 5/22/2014.<br>\n",
      "Adapted (with significant improvements and fixes) from the CVX example of the same name, by Joelle Skaf, 4/24/2008.\n",
      "\n",
      "Topic References:\n",
      "\n",
      "* Section 7.1.1, Boyd & Vandenberghe \"Convex Optimization\" \n",
      "\n",
      "## Introduction\n",
      "\n",
      "Suppose $y \\in \\mathbf{\\mbox{R}}^n$ is a Gaussian random variable with zero mean and\n",
      "covariance matrix $R = \\mathbf{\\mbox{E}}[yy^T]$, with sparse inverse $S = R^{-1}$\n",
      "($S_{ij} = 0$ means that $y_i$ and $y_j$ are conditionally independent).\n",
      "We want to estimate the covariance matrix $R$ based on $N$ independent\n",
      "samples $y_1,\\dots,y_N$ drawn from the distribution, and using prior knowledge\n",
      "that $S$ is sparse\n",
      "\n",
      "A good heuristic for estimating $R$ is to solve the problem\n",
      "  $$\\begin{array}{ll}\n",
      "    \\mbox{minimize}   & \\log \\det(S) - \\mbox{tr}(SY) \\\\\n",
      "    \\mbox{subject to} & \\sum_{i=1}^n \\sum_{j=1}^n |S_{ij}| \\le \\alpha \\\\\n",
      "                      & S \\succeq 0,\n",
      "    \\end{array}$$\n",
      "where $Y$ is the sample covariance of $y_1,\\dots,y_N$, and $\\alpha$ is a sparsity\n",
      "parameter to be chosen or tuned.\n",
      "\n",
      "## Generate problem data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cvxpy as cvx\n",
      "import numpy as np\n",
      "import scipy as scipy\n",
      "\n",
      "# Fix random number generator so we can repeat the experiment.\n",
      "np.random.seed(0)\n",
      "\n",
      "# Dimension of matrix.\n",
      "n = 10\n",
      "\n",
      "# Number of samples, y_i\n",
      "N = 1000\n",
      "\n",
      "# Create sparse, symmetric PSD matrix S\n",
      "A = np.mat(np.random.randn(n, n))  # Unit normal gaussian distribution.\n",
      "A[scipy.sparse.rand(n, n, 0.85).todense().nonzero()] = 0  # Sparsen the matrix.\n",
      "Strue = A*A.T + 0.05 * np.matrix(np.eye(n))  # Force strict pos. def.\n",
      "\n",
      "# Create the covariance matrix associated with S.\n",
      "R = np.linalg.inv(Strue)\n",
      "\n",
      "# Create samples y_i from the distribution with covariance R. \n",
      "y_sample = scipy.linalg.sqrtm(R) * np.matrix(np.random.randn(n, N))\n",
      "\n",
      "# Calculate the sample covariance matrix.\n",
      "Y = np.cov(y_sample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Solve for several $\\alpha$ values"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The alpha values for each attempt at generating a sparse inverse cov. matrix.\n",
      "alphas = [10, 2, 1]\n",
      "\n",
      "# Empty list of result matrixes S\n",
      "Ss = []\n",
      "\n",
      "# Solve the optimization problem for each value of alpha.\n",
      "for alpha in alphas:\n",
      "    # Create a variable that is constrained to the positive semidefinite cone.\n",
      "    S = cvx.Variable(shape=(n,n), PSD=True)\n",
      "    \n",
      "    # Form the logdet(S) - tr(SY) objective. Note the use of a set\n",
      "    # comprehension to form a set of the diagonal elements of S*Y, and the\n",
      "    # native sum function, which is compatible with cvxpy, to compute the trace.\n",
      "    # TODO: If a cvxpy trace operator becomes available, use it!\n",
      "    obj = cvx.Maximize(cvx.log_det(S) - sum([(S*Y)[i, i] for i in range(n)]))\n",
      "    \n",
      "    # Set constraint.\n",
      "    constraints = [cvx.sum(cvx.abs(S)) <= alpha]\n",
      "    \n",
      "    # Form and solve optimization problem\n",
      "    prob = cvx.Problem(obj, constraints)\n",
      "    prob.solve()\n",
      "    if prob.status != cvx.OPTIMAL:\n",
      "        raise Exception('CVXPY Error')\n",
      "\n",
      "    # If the covariance matrix R is desired, here is how it to create it.\n",
      "    R_hat = np.linalg.inv(S.value)\n",
      "    \n",
      "    # Threshold S element values to enforce exact zeros:\n",
      "    S = S.value\n",
      "    S[abs(S) <= 1e-4] = 0\n",
      "\n",
      "    # Store this S in the list of results for later plotting.\n",
      "    Ss += [S]\n",
      "\n",
      "    print 'Completed optimization parameterized by alpha =', alpha"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Completed optimization parameterized by alpha = 10\n",
        "Completed optimization parameterized by alpha ="
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 2\n",
        "Completed optimization parameterized by alpha ="
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Result plots"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "\n",
      "# Show plot inline in ipython.\n",
      "%matplotlib inline\n",
      "\n",
      "# Plot properties.\n",
      "plt.rc('text', usetex=True)\n",
      "plt.rc('font', family='serif')\n",
      "\n",
      "# Create figure.\n",
      "plt.figure()\n",
      "plt.figure(figsize=(12, 12))\n",
      "\n",
      "# Plot sparsity pattern for the true covariance matrix.\n",
      "plt.subplot(2, 2, 1)\n",
      "plt.spy(Strue)\n",
      "plt.title('Inverse of true covariance matrix', fontsize=16)\n",
      "\n",
      "# Plot sparsity pattern for each result, corresponding to a specific alpha.\n",
      "for i in range(len(alphas)):\n",
      "    plt.subplot(2, 2, 2+i)\n",
      "    plt.spy(Ss[i])\n",
      "    plt.title('Estimated inv. cov matrix, $\\\\alpha$={}'.format(alphas[i]), fontsize=16)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<matplotlib.figure.Figure at 0x10ba7bc50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAArQAAALMCAYAAAAVTyF2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3c+WG+d5J+Bfz+E6oqUbENuT/TQjr93HpOcCQtq5gaHi\nCwgVezbhZk7I0QWMKOYCYsnO3pKVA20TWdReFOULiCjaewmzqCoDDQINoBtA9Vv1POf0ARp/qr6v\n8NWLXxUKqAQAAAAAAAAAAAAAAAAAAAAAAAAYsNtJniX5vv17keR3vbboarie5OMkHyX5LMm9JY+5\nddAW1XYzybdJ/kffDenZtmPGcmOdt3O2hj9LU7Pmb/t+y2kuG6d9jMXL1tj7mS2Hh2sea13bnTHU\nueMkHyb5YM39H7WPuXGgdpFmoX+X5M2e23FVfJXkn9vrnyX5zyWP+SjJawdrUW3HSb5MrYK1D9uO\nGcuNTX2bpoYvepgm0P3VFtNaNk77GIu7qrHfZ1bPV7Gu7c6Q69ziTsBfL3nM9TTr40/a/++1z5EX\nDuTDNC/Omz234yq4mWZZdIPxRpYvl23fJMCYYV++yvJAmzSfvL25xbSuyjjdVTs2CbTszlUZP/v2\nfZL/t+T2x3l1J9j3aT4xuPL+W98NYKdeby+/bS+/TvLHhcc8bi+PDtEgBsGY4ZAeJbnTXv8gzV6j\nTVyVcXpV2sF2vG7NHtnfL9z2+yR/30NbtjbEQPtOmkD3fZrjYT5ur3+Ws8eCPM5s13u3Bfwos+Nx\nfzL32Eft8z9K8t7C7d00brTz/iqzrZnbc8/7oL1+ssF0z/Mos+Nj54+tejj3/5Ms/zjhYWbHCHXH\nyNzK7KO9Zf2Yv6/bcv1DVh/fdtE+fdb+LT5nVX83ff2uz/W1O0bvztx0znsN32un1Y2lzrppbjoG\nN+3/uuX5aG5+dzJ7fbqPz7pjFV/k1S3tdX25yJiZX27d6zC/nH/dPqdr57LDYhinH6Z5U522//8i\nyRdz96+qqavG6aqxeOj1pbNufX7UzvfLtq/rPF7Sv328B55nn/V7cVrq3P4ct5ffLNz+p7n72LNl\nhxzcydkV47U0A/jFwnM/yKsv3mc5ezzM4i74D9IM9sV5PUwzeJ+lKUZJM/C7EHijfVw37XXTXebD\nnP3i22IBubUwj2XuZfnHKuf1Y9lz7ufVQHvRPnXz6ZbRP8zdd15/N3n9ui/IdboCNe+8vnfLdL7A\nbjPNdWPwvP5vszy7+X2ZZl04mZv/r9O8dt1rNv8Gu0lfLjJmli23bv7/q/3/nSz/2Itx+SqzL4V1\nYeBvVzz2WVbX1FXjdNlY7GN9Wbc+dyHvzbl5bXLIwXn928V74Hn2Xb8Xja3OXW+nd6e9XPaltetr\n/pZZdsjB7Zx9/TpdxhrDoRi9WxZo7+bVF+beksd1K0O3pXY9Zwf9cc4OzPlpv7nwf7dH9CSzFXJ+\n2kkzIN9cM91V3yjsnjNf6LsVpztguxuQ5xWEt7N8cJ7Xj2XPeSdnC8Jl+tQ952aaAve32ay/616/\npCma80V18fWbv21Z3xePS952mueNwU36v+nyvL3k8R+mOTZx/nVbfIPcpC8XGTPLllvSvD7ft89d\nt7HDOMwfQ3sjTThZFWhX1dRk9ThdNhYPvb6sey+5nldDRlff1gXaZf3b1XvgeQ5RvxeNqc4dpwn4\n8+3pvqB1o21HN/9Vf19muYsE2jcv1o3DudZ3A/bs87nr3y65/2mS50l+leS3aQb0/BbkzfbyF0l+\n3l5/PU0BPs7Z41M/nptm5/dpBkN3+7+2z7l7znRvpDn2dVHXlucL7U+arbZ/W/Kci1jWj02ct6w2\n7dPnSd5or99duG++TbfTvF7rXr+0bXmtnd7tJG+1ty/bct2079tM87wxuEn/t1meSVMA571M8ucV\nj02268sq24yZn6VZB7Y5NpLx+Dpnf07oOM0b92/b/+dr6udpQsYfLzG/Q60v695LXra3fTX3nFXr\n+DYu+x54nkPU71XGUOc+TPJ/crZfz9O0//U0hxYm5x+asrhH/jzda/XGkvumudx6dhBDPIZ23vP1\nD8mjNCvmjTSDfj4Ydl+yeifJ/2z/3kry10n+fWE6iytY0gzqnyV5v53+ozRbyT/YYrqLbdm1k4X/\nl/VjmcVBv82yWnzOtvdN566f9/olTQH7Q5qi8jDn7+3YtO/bTPO8MbhJ/7dZnsusK2jb9KVz0TGT\nNMdjdR+7/u8tnsd4PM4sNNzN2eP35mvqcWY1dZmbK24/zz7Wl5tZvz5397+ce9789Yu67HvgeQ5R\nvzc1tDrXbcgtLo+XSX6Us4eufH7O3x+3mOeqsXL9nPuulKEH2k10WzkfZrYF1ukGzd8s3H43r37s\nu1h8bib5TZotz1+kWcHfT/LTzFaMTabb6dr2w7nbukK/+K3E8yyu+E8W/l9WRJcVi5s5W5i2WVbr\nnnMnm/f3vNfvOM0W8sMk/5Jm5V629dnZ5A1k22me57z+X2R5bmvTvlxkzKzyqJ3n+2ne3BffNOBp\nZnsn/y6zPYvn1dTk1XH6/o7bddH15f2sX5+7+//7wvwO4bwaep591+9dqVjnbubs3vXOUZo9zV/M\nPe7FOX/bfhFtfn3q3MrsFyCutCEE2m53/rKf2vjBhre9n2ZgLL5oT9OsfI8yG5A3k/wyr34ctGz+\nt3L2eK/nSf5jy+l2vk5TzOd/PuPv2zZ3H0mctyw63Ur50zQfVXy1cP+y53ZbZ91Av5mmiB1lFq4u\n0qenaVbaX6VZVt03UafZrL+dVa9ftzy6j5eO0xSX5Gyh7Szr++sL9207zfPG4Hn933Z5rnrtlxXu\nbftykTGzuNyS5s37JM2byi/S7MX4MH60e+xW7c37v2nGy3zQWFVTk9XjdNlYPPT6sm59/lOaenc3\ns/XhH9vLZTVk3rL+dS76Htj9esJ59l2/lxlKnXuc809Y0O1Nnne9ffx80P28beuqv5+tmH6yfJk9\nymzvcNIcBvJVdr9xyIL5s158l2bl+yhNsXvW3vafaVa0e0tumze/hbxM9/Mc3el1u4O0u58J+S6v\nnmb2pH3sw8x+tmTxIOxV0z3P/PTmPza5O9eWZzl70PyiDzL7GZG/WtOPTncqxo/aNnQH+HdfArho\nn15r29M9Z/GA9FX9nXfe6/ewne6zzL4F+1Hb7sWxsdj3+WX6ZWb9XDfNbcbguv5vsjzn29lN+73M\n1o3fpRmPH7f/f5PZsjyvL4tfhNt0zCwutztz7fmybd/ddr7deN30J4IYjvkx1NXwxVPffpfZ2Nik\npi6O02Vjcdn6ue/1pbNufX6vnfbvMvt2/XdZ/jOM2aJ/274HPszqLxTN23f9njekOvdee9t579P3\n2r/uPeUkzbJ6L6sPsznPncx+6aJb57rlMe8kzbL5KLNlBQAAr3gYG/IAABTVHdcLAAAlLZ7JDAAA\nAAAAgLV8xAJQh5pNKX2e+vZOmt99O86rP2A8VN1Pbfwwze8PjsXtNL/t927fDTmQ7sw386cnHLIx\nrstjNMbXWc0eh7HV7GSA63NfJ1boTkn4SXs5hrMF3Urzw9pP0gygxd8BHLLp+ocMyi/TnM3oeoY/\ntk/S/Lj9J+3l0Ps7Vmq2mj1kY6rZyUDrdl+B9ueZnc7weZqtwaE7zqyfz3O40xr27SSzN8ExuJvZ\nj4S/m9k56YfsUXt5nHH0d4zUbDV7qMZYs5MB1u2+Au31nD2t3rJTsA3Nk8x269/M9udYrmrVaS2H\n6q004/kk4zgG7Wma01x2ZwpimNRsNXuoxlazk4HW7b4CbbL8vMhjcDPNeZq/6LshBzC2Lf3Of2W2\nxXvnvAcOwPU0p0+8l+bN/0a/zWGP1OzhU7OHX7OTgdbtvr4U9jKzrcAfpDmf8ljcSvKrvhtxIMft\n3xtpXu+TDOSjjXN8k2bLN2nG+Y/SHJs1VPeSPE7y5zT9vZvxfJFkTNTscVCzh1+zk4HW7b720P46\ns+ORbiT5uKd2HNrbmQ2aMXzB4Lft3zTJaxnHFw1+k9nYvp7kP3psy6H8ub38JE1xZHjUbDV7qMZY\nsxN1e6fupSkQ99Y9cCBupzlW5Vl7+ZN+m8Me3UvzsdU/992QA7mfpr9jWZfHSs1mqMZWsxN1GwAA\nAAAAAAAAAAAAAAAAAOCCdnHmlzH8Th0wUD/+8Y/z6aefjuksWGo2UN0rNbuvM4XtxHRary4/ePAg\nDx486LsZB1W5z0dH/eScPsf22PrcV3+B4Rhjzb5q+jpTGAAA7IRACwBAaQLtgZ2envbdhIMbY58B\ngMMp/aWwisfQUsvYjidNxtfntr9jOghN4YQdG2PN7tkrnbaHFgCA0gRaAABKE2gBAChNoAUAoLRN\nTqxwJ8nLJMdJnuy3OQBckpoNjM66PbQ328tP2suTPbYFgMtRs4FRWhdof57k2/b68yS399scAC5B\nzQZGaV2gvZ7kxdz/b+yxLQBcjpoNjNImXwob5S/2AhSlZgOjs+5LYS+TvN5e/0GSb/bbHID9mkwm\nmUwmfTdjX9RsYJTWbcmfJHkrzTdl7yf5OMkXC49x6lsGa2yngU3G1+eBnfr2StdsGKox1uyebX3q\n26ft5a00W/6LhRGAq0PNBkZpF7HeHloGa2x7K5Px9Xlge2g3oXDCjo2xZvds6z20AABwpQm0AACU\nJtACAFCaQAsAQGkCLQAApQm0AACUJtACAFCaQAsAQGkCLQAApQm0AACUJtACAFDaLk4APO35HOwH\n1+c5mzmskZ4je4zG9EIrYEB1r9Rse2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA\n0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUA\noDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gB\nACjtWt8NuIzpdNrLfI+OjnqZb9Jfn8fK8h6+PtdnYLfGmAv6ctWWtT20AACUJtACAFCaQAsAQGkC\nLQAApQm0AACUJtACAFCaQAsAQGkCLQAApQm0AACUJtACAFDaJqe+vdde/jDJL/fYFgAuT80GRmfd\nHtpbSX6f5EmS4/Z/AK4mNRsYpXWB9jjJ7fb68/Z/AK4mNRsYpXWHHDyZu34zyb/usS0AXI6aDYzS\npl8Ku5nkD0m+2GNbANgNNRsYlU2+FJY0x2H9atWdDx48+Mv109PTnJ6eXqpRAPsymUwymUz6bsa+\nnVuzAarYtGYfbTCtt5O8316/leSThfun0+l0q8ZVd3S0yWLbj7Eta9i3dn3ub6XevbU1+7DNgcPp\n6z2yz1zQl56X9SsLfN0hB7eTPEzyLMmLKIQAV5maDYzSLjYp7KE9oLEta9i3Ae6hXUcRYbDsoT2c\nantoAQDgShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAF\nAKC0a303oKK+zl+c9He+6D77DMB2xvg+1ee8x/ge2efrvIw9tAAAlCbQAgBQmkALAEBpAi0AAKUJ\ntAAAlCbQAgBQmkALAEBpAi0AAKUJtAAAlCbQAgBQmkALAEBpAi0AAKUJtAAAlCbQAgBQmkALAEBp\nAi0AAKUJtAAAlCbQAgBQmkALAEBpAi0AAKUJtAAAlCbQAgBQmkALAEBpAi0AAKUJtAAAlCbQAgBQ\nmkALAEBpAi0AAKUJtAAAlHat7wawnel02st8j46Oeplv0l+fAS5rjDW7T2Nc3t4jG/bQAgBQmkAL\nAEBpAi0AAKUJtAAAlCbQAgBQmkALAEBpAi0AAKUJtAAAlCbQAgBQ2jaB9v7eWgHArqnZwGhsGmhv\nJ/npPhsCwM6o2cCobBponSgYoA41GxiVTQLtSZJP9t0QAHZCzQZGZ5NA+/reWwHArqjZwOisC7S2\n9AHqULOBUbq25v7j9u+NNFv9J0meLj7owYMHf7l+enqa09PTnTUQYJcmk0kmk0nfzdiXjWo2wNAc\nbfi4e0neSfKzJF8s3DedTn3/YOiOjjYdKrtnfLFP7djub4Dvx7k1+/DNGa++6lefNbtPY1zeY+xz\nltTsXbRGoB2BMa6sjMNAA+15rFAHNNKw0ZsxLu8x9jlLarYzhQEAUJpACwBAaQItAAClCbQAAJQm\n0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJR2tINpTKfT6Q4mA8sdHe1i\nmF7MGMd2n8u7R2Pq9PgGNTA0r9Rse2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA\n0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUA\noDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gB\nAChNoAUAoLRru5jI0dHRLiZTxnQ67bsJo9Ln8u5rbI9xjPXV57HVrzFSQw5rjOuU5d0/e2gBAChN\noAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChtkzOF3UxyI8nrSZ7stzkA\nXJKaDYzOJntof5nkt0muJznZb3MAuCQ1GxiddScAvptmS//dcx4zupPO93nOZg7LedgPp68+t/0d\nysnQ1ewlrE+H1Vef+2R5H9wrnV63h/atJG+k2cq/v48WAbAzajYwSpsccvBfSZ621+/ssS0AXJ6a\nDYzOui+FfZPk6/b6yyQ/SnNsFkBJk8kkk8mk72bsi5oNjNK6Ay9upDkm6900H199leTfFh7jeCwG\ny/Fvh+MY2p1Qs5ewPh3WGI/ptLwPbutjaL9Os5V/J81PwCwWRgCuDjUbGKVdxHpb+wyWvSuHYw/t\nwYyugFmfDmuMewwt74Pbeg8tAABcaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIA\nUJpACwBAaQItAAClCbQAAJR2bRcT6fkc7LA3YxzbY+zz2IztNbY+jccYl/cY+7yMPbQAAJQm0AIA\nUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQA\nAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQIt\nAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUNrRDqYx3cE0SplOR9flXh0d7WKY1tLn\nGOtreffV57a/YxpkCtgBjfH9Qs0+LDW7YQ8tAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm\n0AIAUJpACwBAaQItAAClXdvgMXeSvExynOTJfpsDwCWp2cDorAu0J0meJ3k69//T1Q8HoEdqNjBK\nmxxy8Ki9PI7CCHDVqdnA6KwLtE+TfJ3kRfsHwNWlZgOjtC7QXk/yLMm9NMdi3dh7iwC4KDUbGKV1\nx9DeS/I4yZ/TfMngbpJ3990ogH2ZTCaZTCZ9N2Nf1GxgUDat2Udr7r+fs8Ww2+qfN92qZQMwnY6u\ny706Olo3TIenzzHW1/Luq89tf4cyyNTsK2aM7xdq9mGp2e3tGzz3fppvzb6e5T8BM7q1dYwFqk+K\n42EpjuWp2VfIGN8v1OzDUrPb23cw7dGtrWMsUH1SHA9LcRw8BeyAxvh+oWYflprdcKYwAABKE2gB\nAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoLRru5hIz+fz\nhb0Z49geY58Ztr7GdGJcj4Wa3T97aAEAKE2gBQCgNIEWAIDSBFoAAEoTaAEAKE2gBQCgNIEWAIDS\nBFoAAEoTaAEAKE2gBQCgNIEWAIDSBFoAAEoTaAEAKE2gBQCgNIEWAIDSBFoAAEoTaAEAKE2gBQCg\nNIEWAIDSBFoAAEoTaAEAKE2gBQCgNIEWAIDSBFoAAEoTaAEAKE2gBQCgtF4D7WQy6XP2vdDncdBn\nGA5jmyEa2rgWaA9Mn8dBn2E4jG2GaGjj2iEHAACUJtACAFDa0Q6mMUny4x1MB6APnyY57bsRBzSJ\nmg3UNbaaDQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAwLC8neRZku/bv2dJPlu47fstp3lryW03k3yb5H9cuKXbW9aO\nbdzPbDk8XPPYPvo3VNu+btWW/TtJPkiznn2W5E6/zaEYNXs1NbsfQ6/ZSXKc5MM0tZsr7tsk3y25\n/WGa4vBXW0zroySvLdx2nOTLHHYAL2vHRXyf5J/XPKaP/g3Vtq9bpWX/OMn/mvv/TprxJdSyLTV7\nNTX7sIZcs2/n7Mbir/ttDpv4KsuLY5K8SPLmFtPatpjuy67asUlxZHeuyvjZh+/zauH/Ns0eNtiG\nmn3+dNTsw7kq42ffvk/y//puxGX8t74b0JNHme01+iDJ9Q2f97i9PNp5i7ZzVdrBdob8unXr0B8W\nbn8tyY0Dt4XhUbPpg9eNK2d+a/+Habbw/3bFY2+n2aP0UWbHAp6k+air2zX/Uft3K8l77fS+T/KT\ndhqP0uyZ6j5u/cPc817L7DixF2mOi5p3Pc2xLB9ldvzY/Ee2q9rReTTX/veW9O9RO98v276u29p/\nvKR/78z171aSj9vrn+VseHmc2UcZ3Twetf+/mJveOl2fPlvSp0dp+vpZzh5XdpF59/G6dR+hfp9m\n2b2TZrzez+qxNf/x0I25dv7nkj4d0kcLbThO067f9dMcClOzz96vZqvZ+1Z+D+1YfJXZFwy6gbWq\nOD7L7OOFG+1ju2Nh7mX5xw+3cnYAJ7PjB79M8/HYSWYr5q/badzPrMB0Pk6z4nS6lWfeqnY8ztkV\n5IOFaXUF4825eW3y8dV5/esKzWtpCsKLhed+kOSbhds+y+bHF32YZhkms9fjH+bumw9Li8XzIvPu\n43Xr5vkwTR+fZdbnZcu+m393vOo7uXghut5O7057uewLENfX/K3ycZpQsumbIHTU7Nn9araaPW9f\nNVugLWJ+a/9Gzt/aX/wSy/3MisnbWT64b+bVAdxtSc9/SebDth3zz18sTh/k7Ap/N2cL2qp2dHvD\n5uc3/9zreXXAdiv+uuK4rH/dtP9h7rZu5Z9vazePbplez9kicp7FPt1MU+z+du6++dexKxrdcZwX\nmfehX7f5x3YH5J9kVsCXLfu0/fi+fe6my3PRcZo3i/n2PMvsMIGTufmv+vsyy3V9XbWewXnUbDV7\n03mr2bup2eUD7bW+G9CDr3P2pymO0wyC37b//z7NypAkn6cZsH+8xPwWvxDzMsmfz3n8z9MMzrtp\nVtS32tvXHTN2s738RTuNJHk9zRvDcTvftP93vl4zzU18Pnf92yX3P03yPMmv0izjt7P8Y7Vluj49\nn5vXG+31uwv3dfNKmuX220vO+1Cv27yP28un5z6q8bM043Kb4wkXfZjk/+Rsv56naf/rSZ60t90+\nZxqLe3eS5nV7L82y+LcLtg06avaMmr2amn3xmj0IY/1S2OPMBuDdNMWj87P27/329kdptmKXubni\n9vOsG0x30xxfcz3NxxmbfJv1ZprBnDQfZfzP9u+tJH+d5N/n7n8597z56xf1fP1D8qht4400RWTT\nkPP6Be+b7mDei/bxup0s/L/NrwH8KbOPKv/3Fs/rdKFgcXm8TPKjnP0Y9PNz/v64ZLq/T1NQu2m/\nc4H2wTw1+9XrF6VmN9TsgRlroH2a2Zbu32W2lXozyW/SbBn+Is0K+H6Sn7b3L64g7++4Xcdptt4e\nJvmXNAPvjSWPW9aObjD/zcJ9d9MUhu7+/74wv0Potho/zGyLdhOr+nRnbjo/nLu968/vdzDvbVz0\ndXuy8P82b1aP2nm+n+YNcbHQrnMzZ/fUdI7S7LX4Yu5xL875my+i3ceDd9O8IXfe3rJtsEjNns3v\nENTsxtBr9qCMJdCu2jL8v2kG1fygvZWzx2M9T/If7fVu8P40zR6o7qOgbvrzP+1xfcltyfKVpntM\n95zuo4/jzPZuzReBZe14mqYoPMpsRbmZ5Jdp3gj+lKbw383seKV/bC9/sKRN85b1r7Psuctue79t\nT/czKN03cc/zNM0K/Ks0r0v3rdRpmj79Jsnfzz3+79vpLn7MtDjv8xz6dVs2vXnLlv3dNK/xv6R5\nE/9TmuUy/xuwjzM7tmqZbs/EvOvt4+eL5udtW1f9/WzusU/S9PHv2vk/zv7ekBg2NVvNVrPP2kfN\nXrRsmXFFdD/Z8V379yKvnkZx/lvYJ2kOEn+Y2c+KLB4k/UFmP7fxV2kGajePL9MU1jtzt/1nZj8V\n083vd+28um+Bf5PZxx0P23Y+y+wbmh+1j5k/yHyxHZ3uZ0NetPNZPJj9vXbav8vsm5rfZfUZQjbt\n370lt82b3+PQ9XPVwenzXmv72vXpHxbun3+tVn1ktDjvVeb7eqjXbX6MfpazH5UuW/Zde75s23e3\nne937WO7eb3X3jb/RYlF99q/LhCcpFlW72X1R7ardF+S+C5nv4DQLUvYhJqtZi+b9ypq9sVrdtpp\ndL900a1z3fIAuDIexk9mAVShZgMs6I4RA+DqU7MBllg8Kw4AV5eaDQAAAAAAAAAAw+X4IIA61GxK\nudbjvO+k+dHi47x69o2h6n4n7odpfjx7LG6n+WHqd/tuyIF0p22cP7f2kI1xXR6jMb7OavY4jK1m\nJwNcn/s6U1h3Pu1P2sttTwFX0a00Z4V5kmYALf6I9ZBN1z9kUH6Z5lSc1zP8sX2S5sxMn7SXQ+/v\nWKnZavaQjalmJwOt230F2p9ndi7u52m2BofuOLN+Ps/hzsndt5PM3gTH4G5mZ7h5N82pIIfuUXt5\nnHH0d4zUbDV7qMZYs5MB1u2+Au31nD0n9BjOH/wks936NzOe04GuOif7UL2VZjyfZBzHoD1Nc472\n7jSXDJOvQd78AAAMw0lEQVSarWYP1dhqdjLQut1XoE2Sox7n3aebSf6Q5Iu+G3IAY9vS7/xXZlu8\nd/psyAFcT3Pu73tp3vxv9Nsc9kjNHj41e/g1Oxlo3e7rS2EvM9sK/EGSb3pqRx9uJflV3404kOP2\n7400r/dJBvLRxjm+SbPlmzTj/Edpjs0aqntJHif5c5r+3s14vkgyJmr2OKjZw6/ZyUDrdl97aH+d\n2fFIN5J83FM7Du3tzAbNGL5g8Nv2b5rktYzjiwa/yWxsX0/yHz225VD+3F5+kqY4Mjxqtpo9VGOs\n2Ym6vVP30hSIe+seOBC30xyr8qy9/Em/zWGP7qX52Oqf+27IgdxP09+xrMtjpWYzVGOr2Ym6DQAA\nAAAAAAAAAAAAAECfdvFD2WP4WQ9goH784x/n008/HdNJA9RsoKxVNbvPM4WVNZ1OL/z3T//0T5d6\nfsW/y/Z5jPp+zcY0tj/99NO+X+7R6Pu1HtvY1l99HmKfV9VsgRYAgNIEWgAAShNoD+z09LTvJhzc\nGPs8Rl5nhmpsY3ts/U30eQh8KewCptPRdblXR0dj+r5Owxg7nHZ8jWmQ9Ta4jGvgslbVbHtoAQAo\nTaAFAKA0gRYAgNIEWgAASru2wWPuJHmZ5DjJk/02B4BLUrOB0Vm3h/Zme/lJe3myx7YAcDlqNjBK\n6wLtz5N8215/nuT2fpsDwCWo2cAorQu015O8mPv/jT22BYDLUbOBUdrkS2Fj+sFxgOrUbGB01gXa\nl0leb6//IMk3+20OAJegZgOjtO5XDn6d5K00XzC4keTjvbcIYI8mk0kmk0nfzdgXNRsYlE1r9iYf\nTd1L8+WCVT8BM7qTczsf+WG1520eFWPscFadF7ywK1uzjWvgslbV7F0U8dFVKEX5sARa9mmAgXYd\ngRYoa1XNdqYwAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gB\nAChNoAUAoLRrfTcAeFV7rupRmU6nfTdhNPpa1n2Na2MLhs8eWgAAShNoAQAoTaAFAKA0gRYAgNIE\nWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0\ngRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAo\nTaAFAKA0gRYAgNIEWgAASru2i4lMp9NdTKaMo6Oj3uY9tmWdjLPPMER9rctqNgyfPbQAAJQm0AIA\nUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUNomp769117+MMkv99gW\nAC5PzQZGZ90e2ltJfp/kSZLj9n8AriY1GxildYH2OMnt9vrz9n8AriY1GxildYccPJm7fjPJv+6x\nLQBcjpoNjNKmXwq7meQPSb7YY1sA2A01GxiVTQPtrSS/2mdDANgZNRsYlU1+5eDtJO+2128l+WTx\nAQ8ePPjL9dPT05yenu6gaQC7N5lMMplM+m7GPqnZwGBsWrOP1tx/O8kHSV4keT3J3ST/vvCY6XQ6\nvUAT6zo6WrfY9mdsyxr2rV2f+1upd0vNXkLNhuFYVbN3sZYrjgc0tmUN+zawQLsJNfuAxrasYd9W\n1WxnCgMAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRa\nAABKu9Z3Ayrq89zcfZ2T3PnIgarUbBg+e2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSB\nFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChN\noAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABK\nE2gBACjtWt8NYDvT6bSX+R4dHfUy36S/PgNclpoNh2EPLQAApQm0AACUJtACAFCaQAsAQGkCLQAA\npQm0AACUJtACAFCaQAsAQGkCLQAApW0TaO/vrRUA7JqaDYzGpoH2dpKf7rMhAOyMmg2MyqaB1omZ\nAepQs4FR2STQniT5ZN8NAWAn1GxgdDYJtK/vvRUA7IqaDYzOukBrSx+gDjUbGKVra+4/bv/eSLPV\nf5Lk6eKDHjx48Jfrp6enOT093VkDAXZpMplkMpn03Yx9UbOBQdm0Zh9tOL17Sd5J8rMkXyzcN51O\nff9g6I6ONh0qu2d8sU/t2O5vgO+Hmj1yajZDtapm72LEK44joDgyVAMNtOdRs0dAzWaoVtVsZwoD\nAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNo\nAQAo7douJnJ0dLSLyZQxnU77bsLB9dnnPsdXX/0e2zoF7JaazdjYQwsAQGkCLQAApQm0AACUJtAC\nAFCaQAsAQGkCLQAApQm0AACUJtACAFCaQAsAQGkCLQAApQm0AACUJtACAFCaQAsAQGkCLQAApQm0\nAACUJtACAFCaQAsAQGkCLQAApQm0AACUJtACAFCaQAsAQGkCLQAApQm0AACUJtACAFCaQAsAQGkC\nLQAApQm0AACUJtACAFCaQAsAQGkCLQAApV3ruwGwznQ67W3eR0dHvcy3zz6PTV+vMQyVmk0f7KEF\nAKA0gRYAgNIEWgAAShNoAQAoTaAFAKA0gRYAgNIEWgAAShNoAQAoTaAFAKC0Tc4UdjPJjSSvJ3my\n3+YAcElqNjA6m+yh/WWS3ya5nuRkv80B4JLUbGB01p30+G6aLf13z3nM6E5g7JzN4+G84MPXvsb9\nvNC7t1HNNr4YKjV7+FbV7HV7aN9K8kaarfz7u28WADukZgOjtMkhB/+V5Gl7/c4e2wLA5anZwOis\n+1LYN0m+bq+/TPKjNMdmAZQ0mUwymUz6bsa+bFSzHzx48Jfrp6enOT09PUDTALa3ac1ed7DJjTTH\nZL2b5uOrr5L828JjRnfgiGNlxsPxWMM3sGNoN6rZxhdDpWYP30WPof06zVb+nTQ/AbNYGAG4OtRs\nYJR2sSkzus0SW2LjYWt/+Aa2h3YT9tAyWGr28F10Dy0AAFxpAi0AAKUJtAAAlCbQAgBQmkALAEBp\nAi0AAKUJtAAAlCbQAgBQmkALAEBpAi0AAKUJtAAAlHZtFxPp6xzGfZ2zmcPq83Ue49geY5+BYRhj\n/eqrz1eNPbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIA\nUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQA\nAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUNq1vhsAvGo6\nnfY276Ojo17m21ef++ovMBxqdv/soQUAoDSBFgCA0gRaAABKE2gBAChNoAUAoDSBFgCA0gRaAABK\nE2gBAChNoAUAoLRNzhR2J8nLJMdJnuy3OQBckpoNjM66QHuS5HmSp3P/P139cAB6pGYDo7TJIQeP\n2svjKIwAV52aDYzOukD7NMnXSV60fwBcXWo2MErrAu31JM+S3EtzLNaNvbcIgItSs4FRWncM7b0k\nj5P8Oc2XDO4meXfxQQ8ePPjL9dPT05yenu6sgQC7NJlMMplM+m7GvqjZwKBsWrOP1tx/P2eLYbfV\nP286nU63atyuHB2ta/5+9NXfserrdU7G+VqPbb1q+9vfINutK12zgd1Ts9vbN3ju/TTfmn09y38C\nRqBlrwTawxrbejWwQJtc4ZoN7J6a3d6+g2kLtOyVQHtYY1uvBhho1xFoYUDU7IYzhQEAUJpACwBA\naQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClXeu7AZfhfOQM\nVV/n5k56Pz83AFtQsxv20AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm\n0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAaQItAACl\nCbQAAJQm0AIAUJpACwBAaQItAAClCbQAAJQm0AIAUJpACwBAab0G2slk0ufse6HP46DPDNUYX+ex\n9Xls/U30eQgE2gPT53HQZ4ZqjK/z2Po8tv4m+jwEDjkAAKA0gRYAgNKOdjCNSZIf72A6AH34NMlp\n3404oEnUbKCusdVsAAAAAAAAAAAAAAAAAAAA4Er5/2UGlvbQb2jGAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c814dd0>"
       ]
      }
     ],
     "prompt_number": 3
    }
   ],
   "metadata": {}
  }
 ]
}